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1 Introduction 

In this paper we concentrate on approximation to the solution of the linear 
evolution equation 

d t u = Lu = (A + B)u, u(0) = M , (1) 

where L, A and B are unbounded operators. 

As numerical method we will apply a two-stage iterative splitting scheme: 

Ui(t) = exp(At)uo + / exp(As)Bm-i ds, (2) 
Jo 

Ui+i(t) = exp(Bt)uo + / exp(Bs)Aui ds, (3) 
Jo 

where i = 1, 3, 5, . . . and Uo(t) = 0. 

The outline of the paper is as follows. The operator-splitting methods are 
introduced in Section [2] and the error analysis of the operator-splitting methods 
are presented. In Section [31 we discuss the semigroup theory, underlying the 
theoretical method. In Section [3J we discuss the error analysis of the iterative 
methods. In Section [SJ we discuss an efficient computation of the iterative split- 
ting method with ^-functions. In Section [S] we introduce the application of our 
methods to existing software tools. Finally we discuss future works in the area 
of iterative splitting methods. 
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2 Iterative splitting method 

The following algorithm is based on the iteration with fixed-splitting discretiza- 
tion step-size r, namely, on the time-interval we solve the following 
sub-problems consecutively for i = 0,2, . . . 2m. (cf. 



^f- = A*(t) + Bet-!®, with a(t n ) = c n (4) 
and c (i n ) = c" , c_i = 0.0, 

^±iM =^ Ci (t) + Bc 4+1 (t), (5) 

with Ci+i(i n ) =c n , 

where c" is the known split approximation at the time- level t — t n . The split 
approximation at the time-level t = t n+1 is defined as c™ +1 = C2 TO +i(t" +1 ). 
(Clearly, the function Ci+\(i) depends on the interval [t n , t n+1 ], too, but, for the 
sake of simplicity, in our notation we omit the dependence on n.) 

In the following we will analyze the convergence and the rate of convergence of 
the method (@|-([5|) for m tends to infinity for the linear operators A, B :X — > X, 
where we assume that these operators and their sum are generators of the Co 
semigroups. We emphasize that these operators are not necessarily bounded, so 
the convergence is examined in a general Banach space setting. 



3 Semi group theory 

In the theoretical part, we deal with systems of operators. Therefore here we 
have to prove that their operators are generators of Co-semigroup. 

This is not trivial and in the following we give the detail to verify that the 
generators based on a graph norm are bounded operators (see ideas in [5])- 



3.1 2x2 Systems 

We deal with two iterative steps and obtain operators of 2 x 2 matrices. 

Let us assume X to be a Banach space and let A, B be generators of a 
Co-semigroups in X, so we define the system: 

Y.< = (*"), (6) 

c -=( A a'b)- y ^ y «' < 7 > 

where C is a 2 x 2 matrix operator. 

We assume that X^ is the domain of A with the graph norm || ■ ||^: 

HsIUHI^II + lbll, (8) 

where g G dom(A). 
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Theorem 1. Assuming A and B are generators of Co semigroup in and X, 
we have a closed operator C as a generator of a Co-semigroup in Ya 

Proof We solve the Cauchy problem: 

jT{t)f = T(t)Cf, T(0)/ = /, / G dom(C). (9) 

and we find 

T(t)/ " (/*exp(Br)Aelp(!l(t-r))drexp(Bt)) (h) ' /:= (h) ' (10) 



for 



6P:=(^), (11) 



And we have: 




Arm f= — ( exp(At) g 

dt [ >J dt V /„' exp(Br)Aexp(A(t - r))dr 5 + exp(Bi) ft 

exp(At)A g 

exp(Bi)A g + /„' exp(Br)A 2 exp(A(t - r))dr g 
+ cxp(Bt)B h 

_ ( cx V (At) \ / A 

~ \f*exp(Br)Aexp(A(t-r))drexp(Bt)J \A B 

for g G dom(A 2 ) and ft G dom(B). 

The same can be shown for T(t)T{s)f = T(t + s)f, f G P. 

Therefore the family {T(t)} t <o is a Co semigroup in Y^, while R(t) = 
f* exp(Br)Aexp(A(t — r))dr is defined on dom(A) and can be bounded in X^ 
for each t > 0. 

Remark 1. We cannot weaken the assumptions to a the closed operator C as 
generator of a semigroup Y = ' ^ S °^ v i° us ' wnnc if we set A = B we 

find 

f exp(At) \ 
^ ' ~ \tAcxp(At) exp(At) ) { ' 

and if X is a Hilbert space and iA is selfadjoint, then R(t) — tAexp(At) cannot 
be extended to a bounded operator in X for each t > 0, unless we restrict it to 
Xa with the graph norm || • \ \a = \ \A ■ \\ + \ \ ■ ||. 
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3.2 N X AT Systems (Generalization) 

Here, we deal with n iterative steps and obtain operators ofnxn matrices. 

Let X be a Banach space and let A, B be generators of a Co-semigroups in 
X, so we define the system: 

/ X-ABA...A \ 
^ABA...B 



Y 



ABA... A 



X A 
X 



(16) 



and 





(A 


o \ 






A 


B ... 




C := 





B A ... 


■ Y ABA...A — * ~Y ABA...A 




\o 


...A B ) 




We assume that H-aba- 


..A is the domain of A with the graph norm 


\H\aba 


■■A 


= \\ABA- 


■Ag\\ + ... + \\Ag\\ + \\g\\, 



(17) 



(18) 



where g £ dom(A) and we assume 



Remark 2. The proof can be done as for a 2 x 2 operators, and we use recursive 
results. 

4 Error analysis 

We present the results of the consistency of our iterative method. We assume for 
the system of operator the generator of a Co semigroup based on their underlying 
graph norms (see the previous Section [3]). 

Theorem 2. Let us consider the abstract Cauchy problem in a Hilbert space X 

d t c(x, t) = Ac(x, t) + Bc(x, t), < t < Tand x e Q 

c(x,0) = c (x) x e SI (19) 
c(x, t) = a(x, t) xedfix [0, T], 

where A, B : £*(X) — > X are given linear operators which are generators of the 
Co-semigroup and cq € X is a given element. We assume A, B are unbounded. 
Further, we assume the estimations of the unbounded operator B with sufficient 
smooth initial conditions (see J7}j): 



\Bexp({A + B)t)u \ \ < k, 



(20) 
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Further we assume the estimation of the partial integration of the unbounded 
operator B (see J^H): 

\\B [ exp(Bs)sds\\ < tC, (21) 
Jo 

Then, we can bound our iterative operator splitting method as : 

\\(Si-e X p{(A + B)T)\\ <Ct\ (22) 

where Si is the approximated solution for the i-th iterative step and C is a con- 
stant that can be chosen uniformly on bounded time intervals. 

Proof. Let us consider the iteration ((I])-© on the sub- interval [t n ,t n+1 ]. 
For the first iterations we have: 

9 t ci(t) =Acx{t), (23) 

and for the second iteration we have: 

d t c 2 (t) =Acx(t)+Bc 2 {t), te(t n ,t n+1 ], (24) 

In general, we have: 

for the odd iterations: i = 2m + 1 for m = 0, 1, 2, . . . 

d t Ci(t) = Aa{t) + Ba-x{t), t e (t n ,t n+1 ], (25) 

where for co(t) = 0. 

for the even iterations: i = 2m for m = 1,2,... 

d t a{t) = Aa-i(t) + Ba(t), t £ (t n ,t n+1 }. (26) 

We have the following solutions for the iterative scheme: 
the solutions for the first two equations are given by the variation of con- 
stants: 

Cl {t) = exp(A(t n+1 - t))c(t n ), t e (t n ,t n+1 ], (27) 

c 2 (t) = exp(Bt)c(t n ) + f£ +1 exp(B(r +1 - s))A Cl (s)ds, t e (t n ,t n+1 }. (28) 

For the recurrence relations with even and odd iterations, we have the solu- 
tions: 

for the odd iterations: i = 2m + 1, 
for m — 0, 1, 2, . . . 

Ci(t) = exp{A(t - t n ))c{t n ) + /* exp(sA)Bc l _i(r+ 1 - s) da, t G (t n , t n+1 ]. 

(29) 

For the even iterations: i = 2m, 
for m — 1, 2, . . . 

d{t) = exp(B(t - t n ))c(t n ) + f* exp(sS)A Cl _i(t" +1 - s) ds, t £ {t n , t n+1 ]. 

(30) 



The consistency is given as: 

For ei we have: 



ci(r) = exp(A)T)c(t n ), (31) 



c(r) = cxp((A + B)r)c(t n ) = cxp(AT)c(t n ) (32) 
+ / cxp(A S )Bexp((t™ +1 - s)(A + fl))c(i n ) ds. 

We obtain: 

|ei|| = ||c-ci|| < || exp((A + B)r)c(t n ) - exp(Ar)c(t n ) \ | (33) 
< Circ(r). 

For e 2 we have: 

c 2 (t) = exp(S)r)c(t") 

+ y cxp(Bs)Aexp((f l+1 -s)i4)c(t n ) ds, (34) 



c(t) = cxp(BT)c(t") + / cxp(Bs)Aexp((t™ +1 - s)A)c(t n ) ds 
+ / exp(Bs)A (35) 
/ exp(Ap)B exp((r +1 - s - p)(A + B))c(t n ) dp ds. 



We obtain: 



| |e 2 || < \\exp((A + B)T)c(t n )-c 2 \\ (36) 
< C 2 T 2 c{t n ). 



For odd and even iterations, the recursive proof is given in the following: 
for the odd iterations: i = 2m + 1 
for m = 0, 1, 2, . . ., 
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for ei we have : 

Ci(r) = cxp(A)r)c(t") (37) 

+ / exp(As) J Bexp((f l+1 - s)B)c(t n ) ds 
,t B+1 -«i 

+ / exp(Asi)i? / exp(s 2 i?)Acxp((T - s\ - s 2 )A)c(t n ) ds 2 ds\ 

+ ...+ 

+ / exp(Asi)B exp(s 2 B)Aexp((T - s x - s 2 )A)uc(t n ) ds 2 dsx + . . . + 

rt^ / -t"+ 1 -E*I 1 lSl 

+ / exp(Asi)B / exp(s 2 S)j4exp((r — si — s 2 )A)c(t n ) ds 2 ds\ . . . ds i: 

c(t) = cxp(Bt) + / exp( J Bs)Aexp((t" +1 - s)A)c(t n ) ds (38) 
+ ...+ 

+ 1 rt n+1 -si 

+ exp(Asi)B exp(s 2 B)Aexp((T - s-l - s 2 )A)c(t n ) ds 2 ds 1 + . . .+ 

+ / exp^s^i? / exp(s 2 B)Aexp((T - si - s 2 )A)c(t n ) ds 2 dsi . . . 

J exp(s 2 B)Aexp((T - Sl - s 2 )(A + B))c{t n )ds t . 

We obtain: 

|N| < \\exp((A + B)T)c(t n )-Ci\\ (39) 
< CT l c(t n ), 
where a — min* =1 {ai} and < a* < 1. 

The same idea can be applied to the even iterative scheme. 

Remark 3. The same idea can be applied to A — VDV B — — v • V, so that one 
operator is less unbounded but we reduce the convergence order: 

||e 1 ||=^||S||r ai ||e ||+O(T 1+ai ) (40) 
and hence 

||e 2 || = K\\B\\\\e Q \\T 1+ai+a2 + 0(r 1+ai+tt ), (41) 
where < a\,a 2 < 1. 

Remark 4- If we assume the consistency of C(r m ) for the initial value ei(t") and 
e 2 (t n ), we can redo the proof and obtain at least a global error of the splitting 
methods of 0{T m - r ). 

In the next section we describe the computation of the integral formulation 
with exp- functions. 
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5 Computation of the iterative splitting method 

In the last few years, the computational effort to compute integral with exp- 
function has increased because of the ^-function, which reduces the integration 
to a product of exp- functions, see [7]. The ideas are also used for exponential 
Runge-Kutta methods, see [BJ. 

As regards computations of the matrix exponential an overview is given in 

For linear operators A,B : T>(F) C X — > X generating a Co semigroup and 
a scalar t 6 H, we define the operator a = tA and b = tB, and the bounded 
operators <po,A — exp (a), 4>q,b — exp(6) and: 



j\ xp ((l~ s ) T A)j^-^d S ), (42) 
J eMtt-s)TB)-j^-ds), (43) 



9k, A 



»k,B 



for k > 1. 

From this definition it is a straightforward matter to prove the recurrence 
relation: 

<Pk.A = jrl + rAfa+u (44) 

<t>k,B = — ] I + TB(f> k+1 . (45) 

We apply equations (|33]l and (T45|) to our iterative schemes (|2^)) and (|3TJ|) and 
obtain: 

Cl (r) = exp(AT)c(t n ) = 0o,Ac(t n ), (46) 

OO 

Pa(r) = 0o,Ac(t n ) + J] B fc ^fc,A, (47) 
fe=i 

where we assume that B is bounded and exp_B = X^fcLo I^ fe - 

For an unbounded operator B we can apply the convolution of integrals, 
exactly with the Laplacian transformation or numerically with integration rules. 

5.1 Exact Computation of the Integrals 

To obtain analytical solutions of the differential equations: 

dtd = Aci (48) 
d t c 2 - Ac x + Bc 2 (49) 



d t c i+1 = Ac i+1 + Bc i+1 



(50) 
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where c(t n ) is the initial condition and A, B are unbounded operators. 

We apply Laplacian transformation of the differential equations respecting 
the unbounded operators, see pQ. 

We use the Laplace transformation for the translation of the ordinary differ- 
ential equation. The transformations for this cases are given in [1]. For that we 
need to define the transformed function u = u(s, t): 

oo 

Ui(s,t) := Ui(x,t) e~ sx dx . (51) 



o 

We obtain the following analytical solution of the first iterative steps with 
the re-transformation: 

ci = ex.p(At)c{t n ), (52) 
c 2 = A(B - A)' 1 exp(At)c(t n ) + A(A - By 1 exp(Bt)c(t n ). (53) 

The solutions of the next steps can be done recursively. 
The Laplacian transformation is given as : 

C! = (Is + Ay 1 ^ (54) 
c 2 = (Is + B)- 1 c 2 + (Is + B)- 1 A£ 1 (55) 
c 3 = (Is + A)- 1 c ()3 + (Is + A)- 1 Ad 2 (56) 



Here we assume the commutation of 
(A- B)- 1 A = A{A- B)- 1 , 

and we can apply the decomposition of partial fraction: 

(Is + A)- l A(Is + B)- 1 = (Is + A)- 1 (B - A)- 1 A 
+A(B-A)- 1 (Is + B)- 1 . (57) 

Here we can derive our solutions: 

c 2 = cxp( J Bt)c(t") (58) 
+A(B - A)- 1 exp(At)c(t n ) + A(A - B)- 1 exp(Bt)c(t n ) . 

We have the following recurrent argument for the Laplace- Transformation: 
for the odd iterations: i = 2m + 1 
for m = 0, 1, 2, . . . 

Ci = (Is - A)- 1 c n + (Is - A)-^ £i_i, (59) 

for the even iterations: i = 2m for m — 1, 2, . . . 

Ci(t) = (Is - B)- X A g(_i + (Is - B)- 1 c n , (60) 
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We develop the next iterative solution C3 as follows: 



c 3 = exp(At)c(t n ) 
+BA(B - A)-hcxp(At)c(t n ) 
+BA(A - B)- X {B - A)- 1 exp(At)c(t n ) 
+BA(A - B)- X (A - B)- 1 exp(Bt)c(t n ). 



(61) 



We apply the iterative steps recursively and obtain for the odd iterative 
scheme the following recurrent argument: 



+BA(B - A)-Hexp(At)c(t n ) 
+ ...+ 

+BA . . . BA(B - A)- 1 ...(B- Ay 1 ^ 2 exp(At)c(t n ) 

+BA . . . BA(B - Ay 1 ...(B- A)~ 1 {A - B)- 1 {B - Ay 1 exp(At)c(t n ) 

+ ... + BA... BA(B - A)- 1 ...{B- A)' 1 (A - B)" 1 exp(Bt)c(t n ). 

Remark 5. The same recurrent argument can be applied to the even iterative 
scheme. Here we have only to apply matrix multiplications and can skip the 
time-consuming integral computations. Only two evaluations for the exponential 
function for A and B are necessary. The main disadvantage of computing the 
iterative scheme exactly is the time-consuming inverse matrices. These can be 
skipped with numerical methods. 

5.2 Numerical Computation of the Integrals 

Here our main contributions are to skip the integral formulation of the expo- 
nential functions and to apply only matrix multiplication of given exponential 
functions. Such operators can be computed at the beginning of the evaluation. 
Evaluation with Trapezoidal rule (two iterative steps). 

We have to evaluate: 

c 2 (t) = exp(Bt)c(t n ) + J t *„" +1 exp(B(t n+1 - s))A Cl {s)ds, t e (t",t" +1 ], (63) 

where ci(i) = cxp(At) exp(Bt)c(t n ). 

We apply the Trapezoidal rule and obtain: 

c 2 (t) = exp(Bt)c(t n ) + \At (B cxp(At) cxp(Bt) + exp(At)B) , (64) 

where ci(i) = cxp(At) exp(Bt)c(t n ) and At = t - t n . 

Evaluation with Simpson rule (three iterative steps). 

We have to evaluate: 



c, = exp{At)c(t n ) 



(62) 



cs(t) = exp{At)c(t n ) + jl 



+1 



exp(A(t n+1 - s))Bc 2 (s)d, 



Is, te(t n ,t n+1 ], (65) 



where ci(t) = exp(A|) exp(Bt) exp(A±)c(t n ). 
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We apply the Simpson rule and obtain: 

C3(t) = cxp(At)c(t n ) + ^At (b cxp(yl^) cxp(Bt) cxp(A^) (66) 
+ AcxpiA^BexpiAj) cxp(B^) exp(A^) + cxp(At)B^j , 
where a(t) = cxp(At) cxp(Bt)c(t n ) and At = t- t n . 

Remark 6. The same result can also be derived by applying BDF3 (Backward 
Differential Formula of Third Order) . 

Evaluation with Bode rule (four iterative steps). 

We have to evaluate: 

c 4 (t) = exp{Bt)c{t n ) + exp(B(t n+1 - s))Ac 3 (s)ds, t e (t n ,t n+1 ], ( 67 ) 

where Cs(t) has to be evaluated with a third order method. 
We apply the Bode rule and obtain: 

c 4 (t) = exp(At)c(t n ) + ^At (jAc 3 (0) + 32 C xp(B^)Ac 3 (j) (68) 

t t 3+ 3+ \ 

+ 12exp(B-)Ac 3 (-) + 32exp(B-)Ac 3 (-) + 7cxp(Bt)Ac 3 (t)\ , 

where c 3 {t) is evaluated with the Simpson rule or a further third order method. 
We have At = t - t n . 

Remark 7. The same result can also be derived by applying the fourth order 
Gauss Runge Kutta method. 

In the next section we describe the numerical results of our methods. 



6 Numerical Examples 

In the next example, we applied our iterative scheme with their underlying nu- 
merical approximations to differential equations. 



6.1 Linear ordinary differential equation 

We deal with the linear ordinary differential equation: 



du(t) = /-Ai A 2 
dt \ Ai A 2 



u, (69) 



with initial condition uq = (1, 1) on the interval [0, T]. 
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The analytical solution is given by: 



, ci -c 2 exp(-(Ai + A 2 )f) , 

u(<)= U c 1 + c 2 cx p(-( a 1 + a ^), ) ' (7D) 



where 



(72) 



2 1 — T 1 

Ci = j— , C 2 = r 2 -. (71) 

1+ Ai 1+ Ai 

We split our linear operator into two operators by setting: 

du(t) _ ( -Ai 0\ /0 A 2 
dt ~ \ Ai 0J U+ \0-X 2 

We choose Ai = 0.25 and A 2 = 0.5 on the interval [0,1]. 

We therefor have the operators: 

-0.25 0\ (0 0.5 \ 

0.25 OJ ' a ~ ^0 -0.5 ^ ' V A) 

For the integration method we use a time-step size of h = 10~ 3 . 
As initialization of our iterative method we use c_i = 

From the examples one can see that the order increases by one per iteration step. 
In Tables [T} [3] we apply the different integration rules to our iterative scheme. 




Fig. 1. Convergence rates from two to six iterations. 
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Iterative 


Number of 


67" 7*1 


GVT2 


StPDS 


^nlift in cr-nfir'Ht inn 






2 


1 


A ^^91 p 009 


A ^910 009 


2 


10 


3.9664e-003 


3.9664e-003 


2 


100 


3 Q904p-004 


3 Q904p-004 


Q 
O 


1 

_L 


7 fi7fifin HO'} 

f .U f UUc-UUO 


7 fi7fifip 00^ 

( .U I UUc-UUO 


3 


10 


6.6383e-005 


6.6383e-005 


3 


100 


6.5139e-007 


6.5139e-007 


4 


1 


4.6126e-004 


4.6126e-004 


4 


10 


4.1883e-007 


4.1883e-007 


4 


100 


5.9520e-009 


5.9521e-009 


5 


1 


4.6828e-005 


4.6828e-005 


5 


10 


1.3954e-009 


1.3953e-009 


5 


100 


5.5352e-009 


5.5351e-009 


6 


1 


1.9096e-006 


1.9096e-006 


6 


10 


5.5527e-009 


5.5528e-009 


6 


100 


5.5355e-009 


5.5356e-009 



Table 1. Numerical results for the first example with the iterative splitting 
method and the second-order Trapezoidal rule. 



Remark 8. Here we see the benefit of higher quadrature rules in combination 
with the iterative operator splitting scheme, see Figure [1] We obtain the best 
result with a fourth order Gauss Runge-Kutta method. Such improved quadra- 
ture rules and the expansion of the integral formulation show that our method 
has considerable computational benefits. 

In the next example we deal with a Schrodinger equations. 

6.2 Radial Schodinger equation (highly nonlinear) 

We consider the radial Schrodinger equation 



where 



Q^ = f(r,E)u(r) (74) 



f(r,E) = 2V(r)-2E+ l -V±}l , (75) 



If we re-label r t and u(r) — > q(t), ([74")) can be viewed as harmonic oscillator 
with a time-dependent spring constant 

k(t,E)=-f(t,E) (76) 

and Hamiltonian 

H= l -p 2 + l -k{t,E)q 2 . (77) 
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Iterative 


Number of 


67" 7*1 


GVT2 


StPDS 


^nlift in cr-nfir'Ht inn 






2 


1 


A ^^91 p 009 


A ^^910 009 

'-t . O O Zi _1_ t! u U ^ 


2 


10 


3.9664e-003 


3.9664e-003 


2 


100 


3 Q904P-004- 


3 Q904-P-004- 


Q 
O 


1 

_L 


7 f»7fifin 00*3 

i .U ( UUc-UUO 


7 fi7fifip 00^ 

( .U I UUc-UUO 


3 


10 


6.6385e-005 


6.6385e-005 


3 


100 


6.5312e-007 


6.5312e-007 


4 


1 


4.6126e-004 


4.6126e-004 


4 


10 


4.1334e-007 


4.1334e-007 


4 


100 


1.7864e-009 


1.7863e-009 


5 


1 


4.6833e-005 


4.6833e-005 


5 


10 


4.0122e-009 


4.0122e-009 


5 


100 


1.3737e-009 


1.3737e-009 


6 


1 


1.9040e-006 


1.9040e-006 


6 


10 


1.4350e-010 


1.4336e-010 


6 


100 


1.3742e-009 


1.3741e-009 



Table 2. Numerical results for the first example with the iterative splitting 
method and third order BDF3. 



We compare different splitting methods with our scheme, which is related to 
a Suzuki's expansion, see [14] . 

In Figure [H we present the comparison between fractional step (FR), Runge- 
Kutta Nystrom (RKN), standard Magnus expansion (M), improved Magnus ex- 
pansion (BM) and Suzuki's expansion (C). 

Here we see the benefit of the iterative operator-splitting method, which can 
be seen as a modified Suzuki's expansion method. 

Remark 9. The benefit of higher quadrature rules in combination with the iter- 
ative operator splitting scheme is related to Suzuki's expansion. We applied our 
scheme and obtain the best result with a fourth order method. Such improve- 
ments based on quadrature rules, expansion of integral formulations show that 
our method has considerable computational benefits. 

7 Conclusions and Discussions 

We have presented an iterative operator-splitting method as competitive method 
to compute split-able differential equations. On the basis of integral formulation 
of the iterative scheme, we analyze the assumptions of the method and its local 
error for unbounded operators. Under weak assumptions we can prove the higher- 
order error estimates. Numerical examples confirm the method's application to 
differential equations and to complicated Schrodinger equations. In the future 
we will focus on the development of improved operator-splitting methods with 
respect to their application in nonlinear differential equations. 
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Table 3. Numerical results for the first example with the iterative splitting 
method and fourth order GauB RK. 
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